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Direct time radio variability induced by 
non-axisymmetric standing accretion shock instability: 

implications for M87 

Hiroki Nagakura^, Rohta Takahashi^ 

ABSTRACT 

We show for the first time the direct time variable radio images in the context 
of shocked accretion flows around a black hole under the general relativistic 
treatment of both hydrodynamics and radiation transfer. Time variability around 
a black hole can be induced by the non-axisymmetric standing accretion shock 
instability (namely Black Hole SASI). Since the spiral arm shock waves generate 
the density and temperature waves at the post shock region, they cause time 
variability in the black hole vicinity. Based on our dynamical simulations, we 
discuss a possibility of detection for the time variable radio images of M87 by the 
future space telescope VS0P2/ASTR0-G satellite. The most luminous part of 
the images is predicted to be near 15 Schwarzschild radii for some snapshots. We 
show that our results are consistent with existing observational data such as time- 
averaged radio spectra, VLBA images and variability timescale for M87. We also 
discuss observations of M87 with millimeter and sub-millimeter interferometers. 

Subject headings: Black Hole physics, hydrodynamics, instabilities, galaxies: 
kinematics and dynamics, galaxies: nuclei 



1. Introduction 

Elucidating the nature of the strong field of gravity around black holes and the physics 
of the accreting plasma plunging into black holes is a major challenge in physics and as- 
trophysics. Thanks to rapid progress in interferometric techniques from the radio to sub- 
millimeter wavelengths, within the next decade the regions in the vicinity of the massive 
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black holes in e.g. our Galactic center, Sgr A*, and the elliptical ga laxy, M87, whose ap- 
parent sizes are few tens of m i cro-arc seconds, wi l l be directly iinaged (iFalcke. Melia &: Agol 



2000 



Broderick fc Loebll2006l : iHuang et al.ll2007l : iMivoshi et al.ll2007l : iDoeleman et al. 



2008 



Broderick fc Loebl l2009[ ) . The VS0P-2/ASTR0-G project, which is a space- VLBI (Very 
Large Baseline Interferometer) mission, will achieve a spatial resolution of 38 /iarc seconds 
at 43 GHz, corresponding to ~ 2QGM/c^ for the black hole in M87 for a black hole mass 
M = 3.2x IO^Mq and distance D = 16Mpc jMacchetto et al.lll997l : [Macri et al.il999l : lTsuboi 
20081) B Here, G is the gravitational radius and c is the speed of light. 



In the central region of M87, the strong jets are observed from radio to X-ray/gamma- 
ray energies, and high resolution images obtained by the VLBA (Very Long Baseline Array) 



show the bright re gion from which jets are launched (jJunor et al.lll999l : iDodson et al.l 12006 



Walker et al 



20081 ) . The observational data in M87 gen erally show time- variability features 



with timescales ranging from days for TeV gamma-r avs (lAharonian et al. 



20081 ) to months for radio, optical and X-ray bands (IHarris et al 



1997 



2006 



Albert et al. 



Perlman et al.ll2003r i. 



These timescales correspond to, or are larger than, the Keplarian time-scale at the innermost 
stable circular orbit (ISCO) of the black hole in M87, from a few days to one week depend- 



ing on the black hole spin. Theoretically , the observational variabil i ty has been exp 



Becker et al. 2008; Broderick fc Loeb 2009 



based on the accretion disk - jets systems (Chakrabarti fc Wiitalll993t iDi Matteo et al. 



Mandal et al 



20081). In particular, 



ained 



20031: 



Chakrabarti &: Wiita 



( I1993I ) have pointed out that variations in emission from acc retion disks could be explained 
by non-axisymmetric spiral shock waves. iMandal et al.l ( 120081 ) also mentioned that the stand- 
ing shock model can explain the timescale of a few months by comparisons with the free fall 
timescale in the sub-Keplerian flows based on the steady state analysis. As in the examples 
given above, the standing shock waves are important for the time variable emissions in the 
vicini ty of black holes. They have attracted m uch attention frorn many researchers in thi s 
field jNakavamalliggj . ll995[]Das et al.l[2003al jbl: Ichakrabarti et allbooi lokuda et al.ll2007h . 



On the other hand, recent theoretical studies show that the accretion shock waves are 
generally unstable against asymmetric perturbations (namely standing accretion shock in- 
stability or SASI), and they might pla y an important role for the quasi-periodic oscillations 
in stel lar mass black hole candidates (iNagakura fc Yamadal 120081). or for supernova explo- 



sions (iBlondin et al. 



20031: 



Scheck et al 



2008 



According to Molteni et al. ( 19991 ) and also Nagakura fc Yamada J2008I . l2009[ l. standing 



Iwakami et al. 



2008 



Marek fc Jankal 120091 ). 



shock waves are deformed by non-axisymmetric perturbations, and eventually form spiral 
arm structures. Although the linear analysis shows that the shock waves are unstable against 



^ S ince the central r egion around the black hole in Sgr A* is smeared out by intersteller scattering at 43 
GHz (|Shen et al.1120051) . the central region will not be spatially resolved. 
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non-axisymmetric perturbations (iNagakura fc Yamadall2008l . |2009| ). these shock waves keep 
rotating and oscillating in the non-linear phase, while the spiral density and temperature 
fluctuations from shocks are successively plunging into a black hole. These shock waves are 
very differ ent from axisymmetri c ones . Acco rding to linear analysis and dynamical simu- 
lations by iNagakura fc Yamadal (120081 . 120091 ). the rotation and oscillation periods depend 
on the location and strength of the shock waves. They may be related with the instability 
mechanism which still remains elusive. The strongest candidate for the mechanism of SASI 
in black hole accretion (namely Black Hole SASI) is the Papaloizou-Pringle type instability, 
which is induced by acoustic-acou stic cycles between a s hock radius and an inner sonic surface 
(or possibly a corotation radius) (iGu fc Foglizzdl2003l : iGu fc Lull2006l : INagakura &: Yamada 
20081 . 12OO9I ). 



The dynamics of Black Hole SASI are affected by the strong gravitational potential. 
This suggests that if the observational evidence of Black Hole SASI can be obtained, it 
can be used as a diagnostic for the dynamics of plasma around black holes. Thus, it is 
necessary to investigate the possibility of detection. In this paper, we perform hydrodynamic 
simulations of the Black Hole SASI and the radiative transfer by taking into account the 
general relativistic effects around a Schwarzschild black hole in order to investigate the 
observational feasibility of the time variable features from M87. We show that our results are 
consistent with the observational features of M87 such as energy spectrum, image morphology 
taken by the VLBA, and variability with timescales of a few months. We also show that the 
time variability induced by spiral waves swallowed by the black hole will be observed by the 
high detection sensitivity and time resolution of the VS0P2/ASTR0-G satellite. 

This paper is organized as follows. In section 2, we explain the methods and models in 
this work. The numerical results and comparisons with M87 are given in section 3. Finally 
we summarize and conclude this paper in section 4. 



2. Methods and Models 

In this paper, we consider non-axisymmetric shocked accretion flows onto a Schwarzschild 
black hole. We briefly explain the numerical procedure and models in this section. 



2.1. Hydrodynamics 



The detailed procedure for the hydrodynamic simulations is described in INagakura fc Yamada 

(I2OO8I ). We employ a F-law equation of state (EOS) p = (F — 1) poc, where F, p and e are the 
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adiabatic index, pressure, and specific internal energy, respectively. In this paper, we use the 
time-dependent simulation data for the standard model (Ml model) in (INagakura &: Yamada 
20081 ). The initial shocked accretion flows are set as steady solutions which are obtained by 
combinations of the adiabatic index and injection parameters such as Bernoulli constant [E) 
and specific angular momentum (A). We set t hese parameters as F = 4/3, E/c^ = 1.004, and 



A X c/GM = 3.43, respectively. As shown inlljul (119851 ): INagakura fc Yamadal (120081 ) . these 
conditions allow multiple sonic points as well as two-possible shocked accretion flow. Since 
it is well known that the inner shock wave is already unstable against axisymmetric pertur- 
bations in the black hole accretion flow, the outer shock solution comes into existence in a 
realistic system. Thus, we only consider the outer shock which is located at r^ii = IG.lGM/c^ 
in this model. 



According to linear analysis (INagakura fc Yamadal 120081 . |2009|), the outer shock wave 
is unstable against non-axisymmetric perturbations in the linear regime. We have com- 
puted this system in the equatorial plane with a multi-dimensional general relativistic hy- 
dro dynamicaJ_cod^_wliich_js_^^ the high-resolution central scheme, originally pro- 
posed by iKurganov fc Tadmorl (l2000[ l. Actually, our simulations have shown that after 
non-axisymmetric perturbations are added, the initial axisymmetric shock wave starts to 
deform and quasi-steady spiral arm structure was realized in fully non-linear regime. The 
dominant mode for deformed shock waves in non-linear regime is a lower mode (m = 1 or 2 
mode) which keeps rotating and oscillating axisymmetrically and semipermanently until the 
outer boundary conditio ns are changed. We refer the rea der for more details of Black Hole 
SASI hydrodynamics to (INagakura fc Yamadall2008l . 120091 ). In this paper, since we are inter- 
ested in only the non-linear regime, we use these data from enough runs on these simulations 
(t ~ 3 X lO^GM/c^) to the end of our simulations (t ~ 3 x W^GM/c^) for further radiation 
transfer calculations. 



2.2. Radiation Transfer 



By using numerical simulation data described above, we perform general relativistic 
radiative transfer calculations by using the ray-tracing method including the usual spe- 
cial and the general relativistic effects such as Doppler boosting, gravitational redshifts 
and light bending. We assume the rnass M, the distance of the black h ole in M87 as 
M = 3.2 X IO^Mq dMacchetto et al.lll997h . D = 16 Mpc jMacri et al.lll999h . respectively. 
Although in this st udy, we assume the Schwarzschild metric, we u se the code of null geodesies 
in Kerr spacetime ( iTakahashil 120041 : iTakahashi &: Watarail 12007). For r and 6 component s 
of the null geodesies, we use semi-analytic expressions given by iRauch fc BlandfordI ( 1l99ll ). 
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and t and (f) components are numerically solved by the Romberg integration. In this pa- 
per, we only consider thermal synchrotron emissivity and absorption in the accretion flow. 
This is because previous studies have sh own that the observed s pectrum at 43 GHz can be 



fltted by a thermal synchrotr on model (jPi Matteo et al.ll2003l : IWang et al.l l2008l : iLi et al. 



20091 : iBroderick fc Loebl 120091 ) . In this study, we assume the emissivity ji, and the absorp- 



tion coefiicient of the thermal synchrotron radiation of the RIAF model (IManmotd 12000 



Yuan et al.l 120031 ) and calculate the intensity Ii, observed at 43 GHz seen by a distant ob- 



server, by integration along the null geodesies as = f e~'^''g^j^dl where dl is the proper 
l ength differential along the geodesies, r,y = / ai^dl is the optical depth and g = I'ohs/t^ 
( iFuerst fc Wul 120041 : IBroderick fc Loebl l2009l ). Here, i^obs is the observed frequency and z/ is 
the frequency at the local rest frame. When calculating this ^^-factor, we use the simulation 
data of the four-velocity and the special and general relativistic effects are included. The 
emissivity and the absorption coefficients are calculated from the density and the tempera- 
ture of hydrodynamic simulations. We assume a vertical distribution of the number density 
and the temperature of the accretion flow obeying the exponential or Gaussian distributions 
with a Full Width at Half Maximum (FWHM) of r sin(7r/4) for simplicity. Since the density 
profile of the hydrodynamic simulations are normalized by mass accretion rate, we choose the 
mass accretion rate so as to reproduce the observed spectrum. Our numerical simulations do 
not solve the electron temperature and the magnetic field which are required to calculate the 
observed images based on the synchr otron radiation. In t his st udy, the electron temperature 
and magnetic fields are treated as in IBroderick fc Loebl (l2009f ) and we assume plasma beta 
/3 = 10 in all regions. 

By using the simulation data and assumptions denoted above, we fit the observed radio 
spectrum of M87. In Figure [H we show the radial profiles of the electron temperature Tg [K] 
for one snapshot (dots) and the time-averaged data {thick dashed line) for our simulation data. 
Note that since electron temperature depends on the azimuthal direction for one snapshot, 
there are several dots at the same radius in this figure. We also show in Fi gure [J the electron 
temperature profile for the RIAF model used in IBroderick &: Loebl (|2009l ) {dotted line). Al- 
though the electron temperature profiles used in this paper become nearly constant near the 
black hole compared to the RIAF model, the absolute values of the electron temperatures are 
roughly similar to those of the RIAF model. In Figure [2], we show the results of the spectral 
fitting of our models {solid lines) to the observed radio spectra of M87 (cross signs). Here, 



we us e the simulation d ata of eight snapshots a nd the observatio nal data in iBiretta et al. 
Jl99lh (1.49 X 10^° Hzl.lSDencer fc Junoil Jl986h (2.18 x 10^° Hz). IPespringre et al.l Jl996f ) 



(8.90 X 10^° Hz) and lBaath et al.Nl992h flO" Hz). At 43 GHz, the ob served spectrum can be 



fitted by assumirig synchrotron emissivity due to therma l electrons (jPi Matteo et al.l 12003 



Wang et al.l l2008l : IBroderick fc Loebl l2009l : iLi et al.l |2009|) . After calculating the images of 
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the intensity observed by the distant observer, we take into account the effects of the 
spatial resolution at 43 GHz of the VS0P- 2/ASTR0-G sate llite by using the u-v coverage 
data created by the simulation tool, ARIS (lAsaki et al.l 120071 ) and calculate the map of the 
spatially-smoothed intensity /vsop-2 = / / 6^^'^''''+^'^ S{u,v)V{u,v)dudv 

where V{u,v) is the visibility obtained by Fourier transformation of ly, and S{u,v) is the 
sampling function calculated from the u-v coverage of the VS0P-2/ASTR0-G observations. 
By performing the calculations for the simulation data, we finally create the time-variable 
images of the central region of M87 observed at 43 GHz by the VS0P-2/ASTR0-G satellite. 



2.3. Summarizing Our Assumptions 



Since we impose many simplifications in this study, here we summarize the important 
assumptions in order for the reader to understand our models. 

Firstly, we neglect the magnetic field in the gas dynamics, although the radiation is 
assumed to be due to synchrotron emission. This assumption is inconsistent, but we ne- 
glect here for the following reasons. In this study, we consider the 'observational feasibil- 
ity' of Black Hole SASI which does not depend strongly on the effects of magnetic field 
to gas dynamics. The details of gas dynamics should be analyzed under the general rel- 
ativistic magneto-hydrodyna mic (namely GRMHD) treatment in black hole space times 
( iMcKinney fc Gammiell2004j ). On the other hand, the effects of magnetic field to Black Hole 
SASI dynamics are unclear at the present time. Addressing these issues are beyond the scope 
of thi s paper. We are interes ted in how the magneto-hydrodynamic turbulence induced by 
MRI (IBalbus &: Hawleylll99ll ) affects SASI, but these are deferred to future work. 



Secondly, we consider the the rmal synchrotron emissivity and the absorption in the ac- 
cretion flow. In the RIAF model (lYuan et al.ll2003l ). synchrotron radiation due to thermal 
and non-thermal electrons is usually assumed. In this study, as shown in Figure [2], the ob- 
served energy spectrum can be fit by the thermal synchrotron at 43 GHz, thus we neglect non- 
thermal components in this study. At 43 GHz, the out f low components can also be negligible 



compared to the disk component (jJunor et al.l Il999l : iDi Matteo et al. 



2006i'). althou gh it should be included in the other energy bands ( iBroderick fc Loeb 



2003 



Dodson et al. 



2009 



Li et al.ll2009l ). Actually, in our calculations, we can fit the observational data near 43 GHz, 
only using the thermal synchrotron components from accretion flows (see. Figure H]). 

Thirdly, our hydrodynamic simulations do not include meridional plane. For this reason, 
we assume the vertical distributions of the matter as the Gaussian or exponential profile for 
now. The consistent treatment of vertical structure for accretion flows is a major challenge 
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for many years, since the multidime nsional singular surface a ppears in the steady ac cretion 
flows onto black hole (but see e.g. ISincell fc Krolik I (jl997| ): ISalmeron et al.l (120071 )). The 
steady solution is necessary as the initial conditions for the investigation SASI in vertical 
dynamics. We will address these issues in order to investigate more detailed time variability 
from Black Hole SASI. 



3. Results and Implications for M87 



Figure [3] shows the results of the observed intensity distributions {top row) and the 
images smeared by the spatial resolution of the VS0P-2/ASTR0-G {middle row) at 43 
GHz. From left to right panels, we show the time evolution of the images at intervals of 
~ one month for the black hole in M87. The scale of the color bar is normalized in order 
that the time averaged images reproduce the observed spectrum as described in the previous 
section. The size of the images is (SOGM/c^)^, which corresponds (3.8 x lO^^cm)^ for M87. 
The spiral waves due to the SASI around the black hole in M87 evolve with spatial scales 
larger than the spatial resolutions of the VS0P-2/ASTR0-G observations. An exponential 
vertical density proflle is assu med for the case sh own in Figure [3l Thus, if shock waves 
exist in M87, as suggested by iBecker et al.l (|2008[ ) and iMandal et al.l ( 120081 ) based on the 
spectral fltting analysis, the time variation features will be detected by the VS0P-2/ASTR0- 
G satellite at 43 GHz. Almost all emission comes from the inside of the shocks, as a result 
of the density and the temperature peaks of the spiral waves. We also show the images for 
the Gaussian proflle of the vertical density structure (Figure H]), and different viewing angle 
i = 45° (Figure |5]) cases. As can be clearly seen in the bottom panels of these flgures, we 
conflrm the observational feasibility of SASI by the VS0P-2/ASTR0-G. 

The most luminous parts of the intensity are predicted to exist at a distance of e.g. 
30 GMj(? from the vicinity of black hole for some snap shots. This is in contrast to the 
calculations using the stationary accretion flow model where the most luminous part exists 
around the blue-shifted part of the accretion flow around the black hole shadow. This 
means that the position of the black hole in M87 can not be determined from the brightness 
distribution observed in a few months. In our simulations, such time-variable signatures 
will be smeared out by taking the average of the images with the time scale longer than 
the SASI timescale, i.e. several years in the case of M87. Before the launch of the VSOP- 
2/ASTRO-G satellite in 2013, the central regions of M87 will continue to be observed with 
existing telescope s at frequencies of 230 an d 345 GHz, with spatial resolutions of 17 and 11 
/xas, respectively (IBroderick &: Loebll2009l ). It is considered th a t at 3 45 GHz, jet emission 
dominates over disk emission, as shown by IBroderick fc Loeb J2009I ) (e.g. Figure 6). In 
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Figure |3l we show that the effective emitting region at 345 GHz is smaller than that emitting 
at 43 GHz. This is because the accreting plasma is optically thin at 345 GHz and then the 
effective emitting region exists around the shadow. Only when the contamination of the 
emissions from the jets will not conceal the emission from the spiral waves due to SASI, it 
has a possibility that the time evolution features due to the SASI in the intensity map or 
the visibility will be detected by the interferometers at 345 GHz. However, our simulation 
shows that the SASI phenomena will be more clearly observed at 43 GHz because of the size 
of the effective emitting region is larger than the case of 345 GHz. 



4. Discussion and Summary 



In this paper, we have investigated the time variable images and time-averaged spec- 
tra at 43 GHz frequency around a Schwarzschild black hole, employing general relativistic 
hydrodynamic and radiation transfer simulations. The Black Hole SASI causes the time 
variability and its effect most obvious near the black hole shadow. We applied our model 
to M87 and all of our results are consistent w ith past obse r yation al data. The variability 



timesc ale around the core region, reported by [Harris et al.l (119971 ) and also iPerlman et al. 



(120031 ) ■ is naturally explained by our Black Hole SASI model. It should be noted that it does 
not mean to exclude the other steady disk models such as RIAF, since the radial profile of 
the time-averaged electron temperature for our model is similar to RIAF. 

As shown in this paper, the dynamical features in the vicinity of M87 black hole are 
expected to be captured by future space telescopes, such as the VS0P-2/ASTR0-G satel- 
lite. In st eady state, low luminous accretion flow r nodels, a broadband energy spectrum i s 
expected dYuan et al.ll2003[ IjoUev fc KuncicI l2007t iBroderick fc Loebll2009[ \u et al.ll2009[ ). 
For more detailed comparisons with these broadband observational data, there are some 
improvements remaining to be made. Firstly, the fr ame dragging effec ts di rectly af f ect th e 
black hole shadow and time variability. According to lWang et al.l ( 120081 ) an dlLi et al.l (120091). 
TeV v ariability from M87 indicates a rapidly spinning black hole (see also iJoUey &: Kuncic 
( 20071)1. The dependence of Kerr parameters for the imaging are currently undertaken 
( iNagakura fc Takahashill2010l ) . Secondly, our curre nt hydrodynamic simu l ations do not treat 
the meridional plane consistently. As discussed in INagakura fc Yamadal ( 120091 ) . the steady 
solutions for axisymmetric shocked accretion flow in the meridional plane should be con- 
structed. Finally, the magnetic fields may play an impo rtant role for SASI d ynamics. These 
effects have begun to attract attentions in recent years (lEndeve et al.ll2008[ ). Since it is well 
known that the re are the stead y shocked accretion solutions with magnetic field in black 
hole spacetime (ITakahashil |2006| ) . we regard that the SASI could naturally take in place in 
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these flows. In this respect, the effects of MRI to SASI should be investigated. The general 
relativistic magneto-hydrodynamic is strongly required in order to shed light on these issues. 
These efforts are also currently under way. Last but not least, the Black Hole SASI model 
could provide an interpretation for time variability observed in other accreting black hole 
systems, such as galactic X-ray binaries and Sgr A*. 

We would like to thank Shoichi Yamada for useful discussions and suggestions. We also 
thank Dr. Gandhi for valuable discussion and proofreading. This research is supported by 
the Grant-in- Aid for Scientific Research Fund of the Ministry of Education, Culture, Sports, 
Science and Technology, Japan [Young Scientists (B) 21740149 (RT)]. 
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Fig. 1. — The radial profiles of the electron temperature Tg [K] for one snapshot (dots) 
and the time- averaged data {thick dashed line) for our sir nulation data. We also s how the 
electron temperature profile for the RIAF model used in iBroderick fc Loebl (120091 ) {dotted 
line). 
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Fig. 2. — The results of the spectral fitting of our models {solid lines) to the observed 
radio fluxes of M87 (cross signs). Here, we use the simulation dat a of eight snapshots and 
the observationa l data in Isiretta et al. ( )l99ll ) (1.49 x 10^° Hz). ISpencer fc Junorl (119861 ) 
(2.18 X 10^° Hzl.lpespringre et al.l h9m (8.90 x 10^° Hz) and lSaath et all Jl992h flO^^ Hz). 
See also IHoI (119991 ) for additional data near 43 GHz. 
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Fig. 3. — The simulated intensity distributions of the spiral waves caused by the SASI 
around a black hole {top row) and the images smeared by the spatial resolution of the 
VS0P-2/ASTR0-G observations at 43 GHz [middle row). From left to right panels, the 
time evolution of the images are shown of time intervals of 2.8 x 10^ [s] ~ one month for 
the black hole in M87. The size of the images is (80GM/c^)^. The bottom panels show the 
images at 345 GHz for the same model as the top panels. The viewing angle is assumed to be 
i = 30°. For hydrostatic equilibrium, an exponential profile of the vertical density structure 
is assumed. 
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Fig. 4. — Same as the top two rows in Figure [3l but with a Gaussian vertical density profile. 
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Fig. 5. — Same as the top two rows in Figure [3l but with a viewing angle i = 45°. 



